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Abstract 

The dynamics of a three level atom in a cascade (or ladder) configuration with both transitions 
coupled to a single structured reservoir of quantized electromagnetic field modes is treated using 
Laplace transform methods applied to the coupled amplitude equations. In this system two photon 
excitation of the reservoir occurs, and both sequences for emitting the two photons are allowed for. 
An integral equation is found to govern the complex amplitudes of interest. It is shown that the 
dynamics of the atomic system is completely determined in terms of reservoir structure functions, 
which are products of the mode density with the coupling constant squared. This dependence 
on reservoir structure functions rather than on the mode density or coupling constants alone, 
shows that it may be possible to extend pseudo-mode theory to treat multiphoton excitation of 
a structured reservoir — pseudo-modes being introduced in one-one correspondence with the poles 
of reservoir structure functions in the complex frequency plane. A general numerical method for 
solving the integral equations based on discretising frequency space, and applicable to different 
structured reservoirs such as high Q cavities and photonic band gap systems, is presented. An 
application to a high Q cavity case with identical Lorentzian reservoir structure functions is made, 
and the non-Markovian decay of the excited state shown. A formal solution to the integral equations 
in terms of right and left eigenfunctions of a non-Hermitean kernel is also given. 

The dynamics of the cascade atom, with the two transitions coupled to two separate structured 
reservoirs of quantized electromagnetic field modes, is treated similarly to the single structured 
reservoir situation. Again the dynamics only depends on reservoir structure functions. As only one 
sequence of emitting the two photons now occurs, the integral equation for the amplitudes can be 
solved analytically. The non-Markovian decay of the excited state is shown for the same high Q 
cavity case of identical Lorentzian reservoir structure functions, and differs from that for the single 
reservoir situation. 

PACS numbers: 42.50.Ct, 03.65.Yz 
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I. INTRODUCTION 



The quantum optical behaviour of atomic systems coupled to a continuum of the quan- 
tized electromagnetic field modes has been studied since the early days of quantum physics. 
The quantum electromagnetic field is a large system, which can be treated as a bath or 
reservoir. In most cases the atom-field coupling constants and the electromagnetic field 
mode density are slowly varying functions of frequency, and the dynamics of the atomic 
system can be treated via Markovian master equations , or equivalent methods such 
as quantum Langevin equations (see e.g. P] for details of these standard methods). These 
techniques are based on quantum electromagnetic field states with no special distinction 
for any particular mode in terms of photon occupation number, such as thermal states or 
broad-band squeezed states. Naturally, if one mode of the electromagnetic field was in a 
special state, such as a large amplitude coherent state (as in the case where the atom is also 
coupled to a laser field), then this special mode and the atomic system would be treated 
as a small quantum system with the remaining modes constituting the reservoir, so that 
Markovian behaviour would still apply for the small system. 

In certain cases however, such as for atoms in high Q cavities or in photonic band gap 
materials, either the coupling constants or the mode density (or both) are no longer slowly 
varying functions and standard Markovian master equation methods are no longer valid (see 
[3|| for a recent review). A number of non-Markovian methods have been formulated, see, for 



example, references in Ref. |4j. These include non-Markovian master equations 



ime-convolutionless projection operator master equation 



J, the 



, Heisenbe rg equa tions of motion 



stochastic wave-function methods for non-Markovian processes 



methods based on the essential states approximation or resolvent operators 



the pseudo-mode approach 
short-time approximations 



J2CJ2 

3H 




2l|, Fano diagonalization [22j (and |4|, |23j), and various 
The last four approaches are easiest to apply, providing 



clear physical insight into the processes involved. 

One such method is that of pseudo-mode theory 0, Q This method was developed 
for the case of a two level atom coupled to a structured electromagnetic field reservoir in 
the vacuum states and was then restricted to single photon excitations of the reservoir. The 
treatment started from the time dependent state vector for the atom-field system, written 
as a linear combination of one photon states with the atom in the ground state and vacuum 
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states with the atom in the excited state. The basis of the method was that the atomic 
dynamics only depended in this case on the behaviour of a single function, the reservoir 
structure function, defined as the product of the mode density and the square of the coupling 
constant. The complex frequencies and residues of the poles of this function in the lower 
half complex frequency plane enabled so called pseudo-modes to be introduced, one for each 
of the finite number of poles. The non-Markovian equation for the complex amplitude of 
the state with the atom excited (and the field in the vacuum state) could then be replaced 
by Markovian equations involving the finite number of pseudo-mode amplitudes together 
with the amplitude for the excited atomic state. The pseudo-modes are originally related 
mathematically to the reservoir structure function, but in some cases their physical origin 
can be explained. For the case of the atom in a high Q cavity, where the coupling constants 
vary rapidly near the cavity resonance frequencies whilst the mode density is slowly varying, 
the pseudo-modes can be interpreted Q in terms of the cavity quasimodes For the 

case of an atom in a photonic band gap system, no pseudo-mode theory is yet available, 
though a treatment in terms of quasimodes can be used to account for the frequency 
dependence of the coupling constants and mode densities. A treatment of superradiance in 
a photonic band gap continuum is based on the idea of replacing the photonic band gap 
system by a pair of degenerate cavity modes coupled to the multi-atom system and with each 
other, one of the modes being also coupled to a Markovian bath. In terms of the treatment 
in jj], such a case would produce a Fano-profile reservoir structure function, with the Fano 
window representing the photonic band gap. The two cavity modes would correspond to two 
pseudo-modes. The problem for photonic band gap situations is that the mode density is 
actually a discontinuous function of the frequency, and thus the reservoir structure function 
would not have a finite number of simple poles, though approximate representations of the 
reservoir structure function in such a form might be found. 

Leaving aside the difficulties associated with the pseudo-mode theory for photonic band 
gap systems, it would be desirable to see if pseudo-mode theory could be extended to cases 
where multiple photon excitation of the structured reservoir is involved, as the original 
treatment 20( only covers single photon excitation. The limitation of current treatments 
to the single photon excitation case has been noted in Ref. [3j, but some work has been 
carried out on cases of multiphoton excitation of the reservoir. Such a multiphoton situation 
would apply if the two level atom was replaced by a three level cascade (or ladder) system, 
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with an initial condition of the atom in the uppermost state and no photons present in the 
electromagnetic field. The two-step decay will generate electromagnetic field states with two 
photons present. Another case of multiphoton excitation occurs for an excited two level atom 
coupled to a defect mode containing one photon, the atom also being coupled to a photonic 
band gap continuum. Beginning with the essential states approximation, a numerical method 
based on replacing the density of modes by a discretised model has been used in this latter 
situation 13] and in the case of the cascade system (3]]]. Similar numerical methods have 
also been used to treat stimulated emission in a photonic crystal 3^|. The cascade system 
case with one transition coupled near-resonantly to the edge of a photonic gap (and the 
other coupled to a flat continuum) has also been treated via the resolvent operator method 
in Ref. [33]. Although the treatment is analytic, this feature results from being able to 
ignore processes in which the two emitted photons are produced in a different sequence — a 
reasonable approximation if the two transition frequencies are very different. However, a 
more general analytical method would be desirable, and therefore we aim to see if pseudo- 
mode theory can be extended to treat the multiphoton excitation case without having to 
make assumptions about the order in which the photons are produced. Whether an extension 
is possible involves first showing that the atomic dynamics only depends on the behavior of 
reservoir structure functions — in a cascade system we would expect there to be more than 
one reservoir structure function, since two coupling constants are present. A next step would 
be to then introduce suitable pseudo-mode amplitudes, based on the poles of the reservoir 
structure functions and to show that Markovian equations apply to these pseudo-mode 
amplitudes. 

The present paper shows (following the approach of Ref. that in the case of a 

three-level atomic system the dynamics is completely controlled by the reservoir structure 
functions, and gives several methods for determining the atomic and field behaviour. These 
methods could be applied both to photonic band-gap and high Q cavity cases, since the 
general equations (j!4H23j) defining the solution only depend on the reservoir structure func- 
tions and not on the specific type of structured reservoir involved. However, as a test, in 
this paper we only apply the methods to a situation involving a single Lorentzian reservoir 
structure function. This situation could apply when both cascade transitions have the same 
frequency and are equally coupled to a single high Q cavity mode. We also are able to 
interpret the results via an equivalent pseudo-mode model. Situations involving photonic 
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mndgaps could be modelled by appropriate reservoir structure functions (see, for example, 



20 



28|) 



Section |H] of this paper sets out the theory of non-Markovian dynamical behaviour for 
the three level cascade system where both transitions are coupled to a single structured 
reservoir. The state amplitudes are determined from solutions to certain integral equations. 
Approaches to solving the dynamical equations, including a numerical determination of the 
excited state probability for a simple case (and its pseudo-mode theory interpretation) is 
presented in section UTTl Section HVl deals with the simpler case of non-Markovian dynamics 
for the cascade system with the two transitions coupled to two separate reservoirs, again 
with numerical results presented for comparison to the single reservoir case. An alternative 
approach to solving the dynamical equations based on non-orthogonal eigenfunction methods 
is set out in Appendices IXlfDl The paper is summarised in section 

II. DYNAMICAL THEORY FOR A SINGLE RESERVOIR 
A. The Hamiltonian 

The model system we will consider has a three level atom, with states denoted |0), |1) 
and |2), coupled to a reservoir of electromagnetic radiation modes (or heatbath) which is 
to be at effectively zero temperature. The bath modes will be described by a density p\, 
frequency and raising and lowering operators a\ and a\. 

The Hamiltonian for the system is given (in the rotating wave approximation) by 



H = h 



ujx\1){1\ + (u 1 + u 2 )\2)(2\ + Y,uxa{a x 



+ E [9xi (4|0)(1| + «a|1)(0|) +gx2 (4li)(2| + a A |2)<l| 

A 



(1) 



where the atomic transition frequencies are io\ for < — > 1 (i.e. between atomic states |0) 
and |1)) and 0J2 for 1 < — > 2 (i.e. between atomic states |1) and |2)); see Fig. [T] The 
coupling of electromagnetic radiation with frequency uj\ to the transition < — > 1 involves 
the frequency dependent coupling constant Likewise, g\ 2 represents the coupling of the 
electromagnetic radiation field to the 1 < — > 2 transition. Real coupling constants will be 
chosen. It is these frequency dependent coupling constants combined with the mode density 
which define the reservoir structure. 
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B. Coupled amplitude equations 



The Schrodinger picture state vector for the three level cascade atom coupled to the 
quantum electromagnetic field may be written as: 

|*(*)> = c 2 e-^ 1+ ^) i |2)|0 A ) + ^c 1A e-^ 1+ ^) t |l)|l A ) 

A 

+ E c 0AAe- 2 ^|0)|2 A )+ £ c OA ,e-^+^|0>|l A V> , (2) 

A A,/i,(A</i) 

where c 2 , c xx , co AA and co Aai are the amplitudes of the various states in the interaction picture. 
The radiation states included are: |0 A ) in which all the bath modes are in the vacuum state; 
|1 A ) in which the mode with frequency ut x has a single excitation, with other modes being in 
the vacuum state; |2 A ) in which the mode with frequency uj x has been raised to the second 
excitation, with other modes being in the vacuum state; and, |1 A 1 A1 ) in which the modes 
with frequencies uj x and uj^ have a single excitation, with other modes being in the vacuum 
state. A convention for an ordered listing of the modes A for the quantum electromagnetic 
field is assumed, so that double sum over A,/i does not lead to a specific state |1 A 1 M ) being 
included twice. 

The initial state vector is assumed to be of the form 

l*(0)> = |2>!0 A >, (3) 

which allows us to explore the non-trivial case of two photons appearing in the reservoir as 
a result of the interaction with the excited atom. 

Substitution of Eq. (J2J) into the time dependent Schrodinger equation gives a closed set 
of coupled equations for the amplitudes in the situation where the initial state is given by 
Eq. (J3J). If we then take Laplace transforms of the the coupled amplitude equations we 
obtain the algebraic equations 

sc 2 {s) - 1 = -iJ2g X2 c lx (s + i(u x -uj 2 )) 

A 

sc lx (s) = -i 9niCo\n( s + ~ ^0) ~~ i 9fiiC0fi\(s + ii^uj^ - ui)) 

fi,(fM>\) /«,0<A) 

-igxiV2c 0X \(s + i{u x - - ig X2 c 2 (s + i(u a - ux)) 
scow(s) = -igxiV2c lx {s + - uj x )) 

scoxn(s) = -ig^xcx^s + i^x-u^)) - ig X xc Xtl {s + - u x )) (A < /i), (4) 
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where the Laplace transforms of the amplitudes are c 2 , cu, coaa and cox/i and the trans- 
form variable is s. These equations for a cascade system coupled to a single structured 
reservoir, are equivalent to those in Ref. j^]. It should be noted that similar equations 
are given in Ref. (33], for the case where the two transitions are coupled to two separate 
reservoirs. This case is treated below in section II VI The two separate reservoirs case leads 
to simpler equations — firstly, because it is assumed that the A, /1 photons are produced in 
just one sequence (for example 1 2) 1 a) — ► 1 1 ) 1 1 a) — ► |0)|1a1 m )), and secondly, because 
states of the form 1 0) 1 2 a) are not present. Whilst these may be a good approximation for 
the single reservoir case when the transition frequencies are very different, the other pro- 
cess ( 1 2) I 0a) — ► |l)|l/Lt) — ► 1 0) 1 1 aI/j,) ) would need to also be included when the transition 
frequencies are similar, such as in a quantum harmonic oscillator or a Rydberg atom. 

Following the approach of Ref. jl9| we change variables to the reduced amplitudes 
62,61*, b \(i,b \x such that 



c 2 = 6 2 

ClA = gX2hx 

coxti = g\2g^ibo\ti (A < n) 

coaa = gx29xJ>oxx ■ (5) 



Thus we have 



sb 2 (s)-l = -i^2gl 2 b ix (s + i(u; x -uj 2 )) (6) 

A 

sb ix (s) = -i 9%fioXvL(.s + ifa - wi)) -* J2 gli a x^x(s + i(uj^- u^) (7) 

n,(n>\) m,(m<a) 

-ig 2 xl V2b 0X \(s + i(uj x - uJi)) - ib 2 {s + i(uj 2 - u x )) 
sb oxx (s) = -iV2b lx (s + i{ux - u x )) (8) 
sb 0Xfi (s) = -ibi\(s + - Up)) - ia Xfi b ltJ ,(s + - u x )) (A < //), (9) 

where 

gxig^ / in s 

axtM = — • (10) 

gx29ii\ 

Analogous equations to (jBHHl) are given below in Eqs. (|4()H42j) for the case of a cascade system 
coupled to two separate reservoirs. 
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C. Reservoir structure functions and integral equation to determine amplitudes 



Next we eliminate bo\n,bo XX by substitution of (JBJ) and Q in (J7J). This gives 



s + i(u>\ - L0 2 ) + 



0?1 



-J S + l[U X + u) ri —u 1 - UJ 2 ) 



+ 







S + l{OJ X + Up - UJi - UJ 2 ) 

Together with Eq. (jUJ) we now have a set of coupled equations for the b 2 (s) and b ix (s). If 
we eliminate b 2 (s) we obtain an equation for the b\ X alone: 



E 



s is + i{u x - w 2 ) + Yl 

V V 



g 2 v i 



s + i(uj\ + uj„ — uji — uj 2 ) 
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X/i 



+s- 



+ 9 



fj, 2 



b Ul (s + - u) 2 )) = -i . 



(12) 



s + i[uj x +oj IJj — u 1 — UJ 2 ) 
It is useful to rewrite this by by dividing by s, using the properties of the Kronecker delta 
function, and substituting from Eqs. (jl(Jj) and to obtain 



s + i{u x - uj 2 ) + Y 



9 2 m 



s + i(uj\ + u v — uj\ — u 2 
1 



b lx (s + i{oj x - lu 2 )) 



s + i(uj x + u;„ — cji — uj 2 ) 



+ g ^\K{s + i{oo,-uj 2 )) = — . (13) 



We note that the terms involving the quantity a Xfl are absent in similar equations in Ref . [3] , 
resulting in their equations for b\ X being easily solvable. As mentioned earlier, the additional 
terms we have result from allowing for photons to be emitted into the single reservoir in 
two different sequences, an effect not present in the two separate reservoir case treated in 



Ref. 



331 ] . In our case, we next convert the sums to integrals, i.e. Ylu — > J dui ^p{ui M ) , where 



piojfj) is the mode density, so that Eq. (fT3*|) can be written in the form of an integral equation 



with 



fM 

A(u x ) 
B(uj x ,u,j) 
C 



A(u x )f(LU X ) + / duj^B{uj x ,uj^)f{uj^ 



b lx (s + i{u x -u 2 )) 

s + i(u x - uj 2 ) + J dtu ri p(uj ri 



C. 



9%i 



PM 9^i 



S + i{LU X + u v — Ul — uo 2 ) 

1 



2 ' 

|_ 9jj2 

gxig iA s + + o; M - - ^ 2 ) s / 



(14) 

(15) 
(16) 

(17) 
(18) 



The quantities f(uj\),A(uj\),B(uj\,uJ l j) and C are of course all functions of the Laplace 
variable s, but for simplicity of notation s is left implicit. The integral equation (pPfj) is a 
Fredholm integral equation of the second kind (see e.g., |34j). Methods of solution for such 
equations include replacing the frequency spaces by grids of points, thereby converting the 
integral equation into matrix equations that could be solved numerically for each value of s. 
We will discuss one such approach in section IIII Al In Appendix O we also discuss a more 
formal method of solving the integral equation, based on the eigenfunctions of the kernel 
B(u>\, ujfj)/A(uj\) and of its adjoint. 

We also find it convenient to write the integral equation in the form 

J(uj\) + J dUpKiux^^ffan) = d(ux), (19) 

where 

d{u x ) = C/AM (20) 

K(u)\, LUfi) = B(u x ,Up)/A(u x ) . (21) 

We note that the coupling constants and mode density appear in the integral equation 
only in the form u pg 2 " . These forms are called reservoir structure functions, as they contain 
all the essential features of the reservoir and its coupling to the atomic system. Specifically, 
the reservoir structure functions that appear in Eqs. (fTfij) and (fTTj) are 

RiM = p(u\)gli, 

R 2 (u x ) = p(u x )gl 2 . (22) 

As the coupling constants are proportional to dipole matrix elements multiplied by the 
square root of the angular frequency, it is clear that the factors g^/g^i and g X \/g X 2 in Eq. 
(JT7J) are independent of the frequencies uj\ and u^. Hence a third reservoir structure function 
involving the factor is not needed. As the dipole matrix elements would essentially cancel 
out, the factor is of order unity. 

In principle, we can solve the integral equation and thus determine the bi\{s + i{uj\ — to 2)). 
Furthermore, the solutions obtain their particular form from just the reservoir structure 
functions, rather than the density of states or coupling constants alone. 

Next we see that in the new notation Eq. © becomes 

1 i 

h(s) = ^gl 2 b ix (s + i(uj x - uj 2 )) 

s sY 
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I i r 

= ~ ~ -J dvxp{u\)gl 2 f{u>d, (23) 

and again the step to obtaining 62(5) just involves using the reservoir structure function 
_R 2 (c<ja). Note again that f(uj\) is a function of the Laplace variable s, so the decay of the 
initial atomic state |2) described by 62(5) is non-exponential in general. 

Finally, we note Eqs. (jHJ) and (0) imply that box^ and &oaa are fully determined once 62, bix 
are known, and a>\^ (Eq. (fTU|) ) introduces no new frequency dependence as it is independent 
of frequency. Thus all the reduced amplitudes b 2 ,bix,box^, and 6 aa can be determined in 
principle from reservoir structure functions. As we will see next, this is sufficient to determine 
the reduced density operator describing the atomic system. 

Note that the non-Markovian methods could be applied both to photonic band-gap and 
high Q cavity cases, since the general equations (jl4H2H|) defining the solution only depend on 
the reservoir structure functions and not on the specific type of structured reservoir involved. 
Markovian results can be obtained under conditions where the reservoir structure functions 
p(oox)g1 l A2 are slowly varying functions of ujx- Certain integrals give a constant term whose 
imaginary parts are the (formally divergent) frequency shifts and whose real parts are the 
decay rates for the states |1) and |2). 

D. Atomic density operator 

The atomic density operator is defined by 

p A = Ti F I*) (*| , (24) 

and it is not difficult to show that 

PA = \b 2 (t)\ 2 12)(2| + (| duxpMgl, \hx(t)\ 2 ) |1><1| 

Thus we see that the atomic operator only depends on the reduced amplitudes 
^2(^)5 &ia(^), box^{t) (A < /i), and the reservoir structure functions. As the former can be 
determined, in principle, from the reservoir structure functions, we see that the behaviour 
of the cascade atom in the structured reservoir is completely determined by the reservoir 
structure functions (for the initial state given in Eq. (J3J)). 
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On the basis of this key result, it would follow that any existing system could be replaced 
by an equivalent system, provided that the reservoir structure functions were the same 
in both models. This is the basis of the treatment of superradiance in a photonic band 



gap continuum 28] , where the photonic band gap system is replaced by a pair of degenerate 
cavity modes coupled to the multi-atom system and with each other, one of the modes being 
also coupled to a Markovian bath. In terms of the treatment in |4] , such a case would produce 
the required Fano-profile reservoir structure function, with the Fano window representing 
the photonic band gap. The two cavity modes would correspond to two pseudo-modes. 

The absence of any coherence terms in the atomic density operator is a consequence of 
the choice of initial state, Eq. (j3J). The choice of a more general initial state (even with no 
photons present) of the form 

|tf(0)> = (c 2 |2) + Cl |l) + Co |0))|0 A ) (26) 

would require the introduction of a more general time dependent state vector than 
that given in Eq. (J2J), to include additional states of the form |0)|0a), |1)|0a) and |0)|1>). The 
amplitudes for these additional states are not coupled to those for the other states included in 
Eq. (J2J. Again, the solutions for these amplitudes just involve reservoir structure functions 
and are analogous to those already discussed in Ref. for the simpler case of a two level 
atom coupled to a structured reservoir. However, as indicated above, the atomic density 
operator would then include coherence terms. 

III. SOLUTIONS FOR THE STATE AMPLITUDES 

The integral equation ()14j) can be solved in different ways. These include: (a) numerical 
methods based on converting the integral equation to a matrix equation, (b) expansions using 
biorthogonal eigenfunctions and (c) expansions such as the Fredholm expansion 

Q- Only 

the first of these methods will be used here, but as the second approach using biorthogonal 
eigenfunctions may be used in later work and has not been used previously in Quantum 
Optics problems, it is included here in Appendices O for completeness. 
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A. Numerical solution of the integral equation: case of Lorentzian reservoir struc- 
ture function 



As an illustration we consider a greatly simplified example of a three-level system coupled 
to a reservoir with structure. The simplest possible case is that for the same Lorentzian 
reservoir structure function associated with both transitions, with all the couplings and 
transition frequencies equal to each other. That is, we choose a single coupling constant g\ 
such that 

gx\ = 9x2 = g\, (27) 

which amounts to both the dipole moment matrix elements for the transitions being equal. 
The atom will also have two equally spaced transitions which are resonant with the reservoir 
structure, 

W\ = uj 2 = • (28) 

We refer only to uq in the following. Thus for the single reservoir structure function, we 
have 

2 m 2 1 

Rl = R 2 = Pxgx = —. ( ^_ Wo)2 + (r/2)2 (29) 
as in Ref. [20|. The parameter Q represents the strength of the coupling and T represents the 
width of the reservoir structure function. This situation would apply for identical cascade 
transitions coupled to a single high Q cavity mode. Cascade transitions in a photonic band 
gap reservoir would be treated via a different choice of the reservoir structure functions. 

Using this expression for the reservoir structure function we can determine the functions 
A(u\), B(u>\, ujfj) and C in Eqs. (jl6H18|) and then the kernel, Eq. (J2T|) . becomes 

K( v = ™ 2 (s + ijuJx - lj ) + r/2)(2s + i(u x + oy- 2u )) 

[UJX, ^ } 2tt ' s (K-^o) 2 + (r/2) 2 )( s + zK + ^-2 Wo ))Q(^-^ )' 1 } 
where Q(u) is a quadratic polynomial such that 

Q(u) = (s + iuj)(s + iu + T/2) + tt 2 . (31) 

For this model we thus have an analytic form for the kernel, but to go further it appears 
that we need to use a numerical method. We could utilize an eigenfunction method, such 
as that of Appendix El but choose a very simple approach to solve Eq. The process 

is simply to represent Eq. ([T§|) matrix equation 

(K + I)f = d , (32) 
13 



where K and I are matrices and f and d are vectors. We then invert (K + I) to solve for 
f. Thus K{ux,uj^ is represented at discrete frequency points, in effect a discrete basis of 
spatial delta-functions, e.g., = K(u\,u^). Similarly, f(u\) and d(u\) are represented 

at discrete frequency points. From the definition in Eq. (|15p. we see that if we introduce 
the function f(u)\,t) (which we denote as f(t) for short) via 

f(u\, t) = exp[-z(^ A - Ld 2 )t] bix(t), (33) 

then f(t) is the function whose Laplace transform is f(s) = f(ou\,s). However, in order 
to obtain the real and imaginary parts of f(t), we will need the separate inverse Laplace 
transforms f r {s), fi{s) which are the Laplace transforms of the real and imaginary parts of 
f(t) = f r (t) + ifi(t). For complex s the latter Laplace transforms cannot be obtained by 
just writing f(s) as the sum of its real and imaginary parts. However, the Laplace transform 
f r (s) of the real part f r {t) is real (and similarly the Laplace transform fi(s) of the imaginary 
part fi(t) is real), if the Laplace transform parameter s is real. Hence, the real and imaginary 
parts of f(s) are equal to the Laplace transforms of the real and imaginary parts of f(t) 
for s on the real axis, so Re/(s) = f r (s), Im/(s) = /j(s) for s real. As f(s) is an analytic 
function of s, the analytic continuation of f r (s) + ifi(s) from the real axis will determine 
f(s) everywhere. 

In this example, if we discretise K(u\, Uy) on an TV x iV grid we define the N x N complex 
matrices K r and Kj from the real and imaginary parts of Eq. (|3T)|) on the real s-axis, and 
then Eq. (|3*2|) becomes 



K r + I 






K r + I 



_cL_ 
d, 



(34) 



for s on the real axis. The formal solution for f r = f r (s) and fj = fj(s) for real s then generates 
the solution for f (s) everywhere. Because of this (and having first identified K r , Kj, d r and 
dj for real s using Eq. (|21j) and Eq. (|20j)). we can now regard Eq. ()34|) as applying for all 
s. This approach could not be used if the the real and imaginary parts of K{uj\,uj^) and 
d(u\) on the real s-axis are not analytic. The matrix inversion step thus involves a matrix 
with 4N 2 elements compared to, say, 0(N A ) elements represented by Eqs. (0 - EJ) in an 
equivalent discretised form. 

Thus we solve for f r and fj in Eq. (|3*ljl . and hence determine the bi\(s + i{uj\ — 0J2)) of 
Eq. (jl5|) . We then find the b2(s) from the scalar product form Eq. (}2"3*|) obtained from Eq. 
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(JUJ) so that 



b 2r (s) = (l + r-f^/s 



6 2r (s) 



(35) 



where r = {p\g\ 2 }- Finally, &2(£) is determined by a numerical inverse Laplace transform. 

Figure 121 shows some results for this numerical matrix approach with the kernel given 
in Eq. (|3(jp. which was derived from the reservoir structure function in Eq. (J2T?)) . The three 
curves show the upper state population for three different sizes of matrix which were used 
to discretize the integral equation. Each case used the same parameters Q=l and T=l, 
where there is a distinct non-Markovian evolution that could not be treated perturbatively 
because of the strong coupling to the reservoir structure. The solid curve in Fig. El shows a 
good result that was obtained with a matrix of size 150 x 150 for this problem. Reducing the 
matrix size to 100x100 (dashed) results in only a slight degradation of the result. However, 
further reduction of the matrix size affects the numerical result quite badly. 

The effect of changing the coupling strength Q is shown in Fig. El The probability 
of finding the atomic system in the highest atomic state is shown. For strong coupling 
(Fig. ISfa)) we see damped oscillations that are a typical manifestation of non-Markovian 
processes. As the coupling is reduced, (Fig. Efb)), the oscillations weaken and then further 
reductions in the coupling strength Q, (Fig. EKc)) result in no oscillations and decay that is 
closer to exponential and on a longer time-scale than the strong coupling cases. Fig. |3fc) 
still shows some visible initial quadratic behaviour because of the relatively high value of 



B. Equivalent pseudo-mode model 

The reservoir structure function given in Eq. ()29)1 is extremely simple and as result we 
can reproduce the results of Fig. 01 i.e. the population |&2(£)| 2 , from the Markovian master 
equation 



n/r. 




(36) 



which is given in the interaction picture with the atom-'field' coupling term 



v = n (a t |o)(i| + a|i)(o| + a t |i)(2| + a|2)(i|) . 



(37) 
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In this master equation we have introduced a single oscillator, or pseudo-mode 20 , which is 
represented by the harmonic oscillator operators a and af. In this approach (see |2(|) pseudo- 
modes are introduced as assumed bosonic entities, rather than via constructing pseudo-mode 
amplitudes. A cascade atom resonantly coupled to a damped high-Q cavity mode, which 
is also coupled to a Markovian bath of vacuum modes, is an example of a physical system 
which has the same master equation as (J3"6*j) . Such a model was considered in our earlier work 
(4], where we showed that multiple excitations of a structured reservoir could be treated for 
reservoir structure functions such as Eq. (|2*U|) . To utilize the present pseudo-mode model we 
solve the master equation (jHfij) with the initial condition of an empty pseudo-mode and the 
atom in the state |2). On tracing out the pseudo-mode, to obtain atomic properties alone, we 
can reproduce the results of the matrix method used with the kernel of Eq. (J3(Jj) . It should 
be emphasised that it does not appear to be easy to find such a simple master equation for 
more complex reservoir structures such as photonic band gap models with branch cuts in 
the reservoir structure function. In such a case the approach outlined in this paper (which 
only depends on the reservoir structure functions) may be useful instead. For the present 
Lorentzian model, the agreement between the matrix method given earlier in this section 
and the master equation (jSfij) is excellent. 



IV. DYNAMICAL THEORY FOR TWO SEPARATE RESERVOIRS 

In this paper we have commented in several places that there are differences in our 
single reservoir treatment from the simpler case of separate reservoirs coupled to the two 
transitions in our model system. In our model, the two photons may be emitted in either 
order, whereas with the distinguishable photons in the two reservoir model, only one order 
of emission is involved. So, with the formalism now complete, it is instructive to look at the 
explicit differences between our model and the simpler two separate reservoirs model of the 
kind considered in Ref. [3^. In this case the Hamiltonian in Eq. is replaced by 



H = h 



1) ( 1 1 + (wi + uj 2 ) 1 2) (2 1 + wa4o a + E U A K 



+ E^i (&llo)(i| + ^|i)(o|) + E^ 2 (ai|i)(2| + a A |2)(i| 



(38) 
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where the bath operators a\ and a A for the first bath now couple only to the 1 < — ► 2 
transition, and the new bath operators bj^ and for the second bath couple only to the 
« — ► 1 transition. For the initial state vector, Eq. ©, the state vector analogous to Eq. 
(J2J) no longer contains a term involving Co AA , and there is no restriction over the double sum 
A,//, since the two types of bath modes are now distinct. We can write 

|*(*)> = c 2 e-^ 1+ ^ i |2)|0 A )|0 At ) + 5:c 1A e-^ 1+ ^|l)|l A )|0 At ) 

A 

+ EcoA M e-^ + ^|0)|l A )|U, (39) 

involving product states of: the atom, and one or zero excitation states of the two baths. 
The equations for the Laplace transforms of the reduced amplitudes, Eqs. dOKJ) are then 
replaced by 



sb 2 (s)-l = -i^2gl 2 b lx (s + i(u x -u 2 )) (40) 

A 

sb~i\(s) = -«E^iW( s + K^p - ^i)) - ih(s + i(uj 2 - u x )) (41) 

sbo\fi{s) = -ib!x{s + i{uJi-uJn)), (42) 



We note that at this point the differences are that, as well as the absence of the bow terms, 
there are no terms involving a AM (as in Eqs. (J7J) and (J9~|)) and there are no restrictions over 
the sum over \i (as in Eq. (JJJ)) These equations are equivalent to those in Ref. 33^ . 

As in the case of both transitions coupled to one single reservoir, the dynamical behaviour 
only depends on reservoir structure functions, and following the same approach as in section 
III Dl it is easy to see that the atomic density operator is also determined from these functions. 

If we now follow the elimination procedure of section III CI we find the same equations 
fll4M18|) for /, A, B and C except that the consequence of no a A/x term being present in Eq. 
(JUJ) is that the quantity B becomes 

BK, w M ) — £?K) = piurf-f (43) 

Crucially B no longer depends on uo\ as previously. Expressions for /, A and C are otherwise 
unchanged. 

The integral equation then simplifies to the easily solvable form 

A(u x )JM + J du,B(u,)JM = C, (44) 
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for which the solution is 

— CI 

fM = l + jd^K^u,) ' %)' (45) 

In this case the equivalent kernel is separable: Kiojx^uj^ = B(uj^)/A(uj x ). 

We can apply our results to the situation analogous to that treated in section IIII Al 
where both reservoirs, though now separate, have identical coupling constants and reservoir 
structure functions, and the two atomic transitions are equally spaced and resonant with 
the reservoir structures. We utilize Eqs. (j27H29j) and, for this simple model, the kernel can 
be easily obtained as: 

K( \-^L. (s + i(u x -u ) + T/2) 

1 A ' »>~ 2tt ' s((u,-u y + (r/2y)Q(uj x -u y l4bj 

This result may be compared to the previous expression in Eq. (|3Uj) for the case of a single 
reservoir. 

The integral / duJ^Kiuj^ u; M ) can be performed by using a contour in the lower-half plane, 
and we obtain 

d^K{u^) = Y (s + r /2)(t+r) + ^- (47) 

We may now find from Eq. (fT7)|) that 

Q 2 

A(uj x ) = s + i{uj x - uj ) + — — — r (48) 

s + 1 1 1 + l\UJ X - UJ ) 

so the solution for f(uj\) can be obtained from Eq. (|43j) . We find that: 

? = Ms + m(s + r) + [s + i(u x - up) + r/2] 

J{ X) " (s + T/2) [s(s + T) + 20?} Q(lu x - u ) ' 1 ' 

A numerical inversion of f{uj x ) can be performed to obtain 62 (t) using the same approach 
as in section UlI Al However, The reservoir structure, Eq. (J2T?|) . is sufficiently simple that 
a solution for 6 2 (t) can be found from Eq. ()49|). We first need to perform the integral in 
Eq. (|2*Hj) which is facilitated by the fact that Eq. (j4*U|) has no poles in the lower-half complex 
plane (for Re(s) > 0), while the factor p{uj x )g\ 2 in Eq. (|23|) has only a single pole in the 
lower-half complex plane if we use the example given in Eq. ()29)) . Then if we perform the 
integral of Eq. (|2l?|) around the single lower-half plane pole we find that 

Hs) = \~ n \{s + Y,2)CJ + T) + 2^Y (50) 
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If we now perform the inverse Laplace transform, we find 

Ut) = f 2 e~^ + (l - oosOJt) + ^e~^ sin(/3t) , (51) 

where (3 2 = 2Q 2 - (T/2) 2 . 

The result for the time evolution of the probability for finding the atom in the highest 
atomic state is seen in Fig. There is clearly a difference from the single reservoir result 
shown in Fig. |3] (the dashed line in Fig. The present situation, where both atomic 
transition frequencies are equal and resonant with the structured reservoir, should highlight 
the difference between the cases of two separate or one single reservoir. In this situation 
both photons emitted should have similar frequencies, and the single reservoir case where 
the first emitted photon cannot be distinguished from the other should give different results 
to the two distinct reservoir case where they can. 

We note that for strongly coupled systems, 2Q 2 > (T/2) 2 , the time evolution in Eq. (JoTj) 
can be re-expressed in the form 

Ut) = 2Q2 _ (r/2)2 sin 2 (/3t/2 + 0)e" rt / 2 , (52) 

where 

T/2 

cos0 = -^. (53) 

What is interesting here are the oscillations which are given by the square of a sine function, 
i.e. the probability oscillates as the fourth power of a sine function which is damped at 
the rate T. In the limit Q ^> T the angle approaches 7r/2 and Eq. (|H2*|) reduces to 
b 2 (t) w cos 2 (nt/V2)e- rt / 2 . 

Conversely, for weakly coupled systems, 2Q 2 < (T/2) 2 , the time evolution in Eq. (joTj) can 
be re-expressed in the form 

b 2 (t) = (r/2 jf m2 sinh 2 ( 7 t/2 + 0e- rt/2 , (54) 



where Y = (T/2) 2 - 2tt 2 and 



r 12 

cosh£ = — (55) 



In the extreme limit of Q T, Eq. (|54|) reduces to the Fermi Golden-rule result: b 2 (t) 
cxp(-2tt 2 t/T). 



19 



V. CONCLUSION 



The dynamical behaviour of a three level atom in a cascade configuration in which both 
transitions are coupled to a single structured reservoir of electromagnetic field modes, and 
initially in the upper state, has been analysed via Laplace transform methods. This situation 
involves a two photon excitation of the reservoir, and our equations take into account the two 
possible sequences in which these two photons are emitted. We have shown that the atomic 
density operator is determined from the solutions of integral equations, in which the prop- 
erties of the structured reservoir only appear via so-called reservoir structure functions, all 
essentially given by the product of the mode density times the square of coupling constants. 
In the cascade system two distinct reservoir structure functions are involved since there are 
two transitions. The dependence of the dynamics solely on reservoir structure functions is 
the necessary condition for treating structured reservoir problems via pseudo-mode theory, 
so our results suggest that it may be possible to extend pseudo-mode theory to problems 
involving more than a single photon excitation of the reservoir. 

This result also shows that any existing system could be replaced by an equivalent system, 
provided that the reservoir structure functions were the same in both models. This is the 



2jj and the 



basis of the treatment of superradiance in a photonic band gap continuum 
general treatment of multiphoton excitation in terms of quasimodes given in our earlier 
work j^]. 

In addition, a similar treatment of the dynamical behaviour of a three level atom in a 
cascade configuration coupled to two separate structured reservoirs of electromagnetic field 
modes, and initially in the upper state, has been carried out. One reservoir is coupled to the 
upper transition, the other to the lower transition. This situation again involves a two photon 
excitation of the reservoir, but now only one possible photon emission sequence is involved. 
In this situation, the equations are simpler, and the integral equation for the amplitudes 
can be solved analytically. Again, the dynamical features only depend on reservoir structure 
functions. 

A numerical method of solving the integral equations based on discretising the frequency 
space has also been obtained, and which can be applied to various structured reservoir 
situations - such as for high Q cavities and photonic band gap systems. Here we have applied 
this method in a numerical test for a high Q cavity situation, where the same Lorentzian 
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reservoir structure function applies for both transitions, showing the non-Markovian decay of 
the excited state. Results for both the single structured reservoir case and the two separate 
reservoirs case have been obtained, showing the different behaviour in the two cases. This 
difference is to be expected, as the two photons emitted should have similar frequencies, and 
only in the two separate reservoirs cases should it be possible to distinguish which order 
the photons were emitted. In this latter case we were able to solve the model problem 
analytically. Finally, a formal solution of the integral equations based on the biorthogonal 
left and right eigenf unctions of the non-Hermitean kernel has been presented for completeness 
in the Appendices. 

Our treatment of the cascade system coupled to a structured reservoir may be compared 
to those of Ref. in the two separate reservoirs case and to Ref. 31 1 in the single reservoir 
case. Both these papers also demonstrate non-Markovian decay of the excited state. Our 
fundamental amplitude equations in sections fll Bl and HVl agree with those of these authors. 
The work in Ref. j^j] differs from our treatment, being based on replacing the structured 
reservoir with discrete modes and then using numerical methods. The work in Ref. [3] is 
analytic. However, a direct comparison of the numerical results is not yet possible with either 
Ref. [31] or Ref. j^J , since both applied their theory to a photonic band gap system whereas 
our present application is for the equally important situation of a high Q cavity. Further 
applications of our theory involving good analytic approximations to the reservoir structure 
functions for photonic band gap systems will, however, enable more detailed comparisons to 
be made. 
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APPENDIX A: INTEGRAL EQUATION KERNEL AND ITS EIGENFUNC- 
TIONS 

The kernel K(u)\, uj^) involved in the integral equation (fT^j) and given by Eq. (}2*Tj) may 
now be used to define an integral operator K. The effect of K on any function is defined 
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by 

(K<i>) ux = J du/JCiux,^^) . (Al) 
The eigenfunctions (f) n (u\) and eigenvalues £ n for the integral operator K then satisfy 

k<f) n = e„0„, (A2) 

or (in full) 



M ■ (A3) 
Note that we are following Ref. [35] in our definition of the eigenvalue of the i nteg ral equation, 
rather than the definition used in many mathematical textbooks (e.g. Ref. 3J|) where l/£ n 
would be the equivalent eigenvalue. 

Similarly to Eq. (jAl|) we can define the adjoint operator via 



(tfty) WA = y d^iT(uv,u; A )0(uv), (A4) 
with eigenfunctions 9 n (u\) so that 

K ] 9 n = CJn ■ (A5) 

It is straightforward to show that has eigenvalues which are complex conjugates of 
those for K (see Appendix IB1 for details). As K will in general be non-Hermitian, the 
eigenfunctions <p n do not satisfy standard orthogonality conditions. Instead the <p n and the 
6 n satisfy so-called biorthogonality conditions 

duj^Kux^miuJx) =5 nm . (A6) 

The normalization result of unity for n = m can be arranged by scaling either the 9 n or m 
by appropriate factors. Although these results are familiar in regard to the mode functions 



for unstable optical systems 
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J), they are not widely used in quantum optics. So, for 
completeness, a derivation of Eq. (|A6|) is presented in Appendix O A formal method of 
determining the eigenfunctions <fr n and 9 n is set out in Appendix El 

APPENDIX B: REPRESENTATION OF THE KERNEL 

We expand (p n in an orthonormal basis u n so that 

<frn(uJ X ) =J2 a m U m(U\) (Bl) 

m 
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with 

J duJxU*(ux)u m (uJx) = Sim ■ (B2) 

Then we can write Eq. (jA3|) as 

£ / = £ n ^a™U m (u; A ) • (B3) 

m m 

Then if we multiply by u*[(u\) and integrate we find 

X] (/ / duJxduj^(uJx)K(uJx,uj^u m (uJ IM ) - 5 lm £^j a™ = 0, (B4) 

which must be true for all values of Z. Equation (jB4|) is a matrix eigenvalue equation with 
the matrix 

Kim = J J duj x duj^(uJx)K(uJx,uj Ll )u rn (uj^ (B5) 
and eigenvalues which satisfy 

\K lm -£5 lm \ = 0. (B6) 

For the operator (see Eq. (jA4|0 the matrix is replaced by its adjoint and clearly its 
eigenvalues are complex conjugates of those for K. 

The explicit form for Ki m is found by substitution of the expressions (|21j) and (|17|) into 
Eq. (jB5|) which yields 

Kim = duj x duj u Ui(uJx) 



1 



/ \ 2 9\l9fJ.2 



+ 



U m {0Jfj) . (B7) 



A(u)x) V ' """^ fl^/d s + i(uj x + w M - wi - u 2 ) s 
The integral over uj^ will involve the reservoir structure functions as defined in Eq. 

The function A(u\) is also obtainable from the reservoir structure functions (see Eq. 

Note that gxig^l gxig^i is independent of frequency in Eq. (|B7|) . 

In summary, the matrix Ki m and hence the eigenf unctions n , 9 n and eigenvalues are all 

obtained from the reservoir structure functions and given functions, such as the basis set u n . 



APPENDIX C: BI-ORTHOGONALITY OF EIGENFUNCTIONS 



To show that the eigenfunctions satisfy a biorthogonality condition we first write from 
Eqs. jUl and (jH El El 

J du^K(uJx,UJ^)(pn(^^) = in4>n{^\) 

duJ ri K(uJ^UJx)e m {uJ^ = £ m m (u\) ■ (CI) 
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After multiplying the first equation by O^ux), the second by 4> n (^x) and then integrating 
over ojx we find that 



dLOxdu ii e* m {uJx)K{ujx 1 0J tl )(j) n {u tl ) = £ n / duj x 0^(ujx)4>n(^x) 
du x du ll (t) n {uJx)K(uJx,uJ ll )e* m {uj^ = £ m J du x d* m {ux)<t>n{^\) ■ (C2) 

After a change of variable in the second equation, the left-hand sides are equal and we then 
conclude that 

(£« - U) J du x e* m (ujx)<f>nM = o, (C3) 

so that the biorthogonality condition 

J dcux9* m M<i>nM = (C4) 

applies unless £ m = £ n . 

APPENDIX D: INTEGRAL EQUATION SOLUTION IN TERMS OF EIGEN- 
FUNCTIONS OF K 

We will assume that the set of eigenfunctions n form a basis for expanding the solution 
f(u>x) (to Eq. (fTHj) ). Likewise we will assume that d(u>x) can be expanded in terms of the 
4> n so that 

70a) = YsJnMuJx) 

n 

d(ux) = X]^«N' ( D1 ) 

n 

Using the biorthogonality of the eigenfunctions (Eq. (|A6p ) the expansion coefficients can be 
found as: 



fn = J dwxQ* n {u\)f{ux) 

d n = [ du x e* n {ux)d{ujx). (D2) 



Substituting from Eq. (|D1|) into Eq. (JT5j) and using the eigenvalue equation ()A3|) we find 
that: 

Y.'fn ( l ) n(uJx)+Yjn d U (u X , U M )0„(^) = ^ rf„0„(w A ) 

ran n 

E(7n + CJn - 4)0nK) = 0. (D3) 
n 
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Using the biorthogonality result for the eigenvalue £ n we see that 



/ n (l + £„)-d„ = 0, (D4) 

so that provided £ n 7^ — 1 

In = T^T, (D5) 

which gives the solutions for the expansion coefficients for f(u\) in terms of known quantities. 
The quantities f n , <f> n (u>\), K(u>\, oj^) and £ n are of course all functions of the Laplace variable 
s, but for simplicity of notation s is left implicit. 
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FIGURES 
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FIG. 1: The three level cascade (or ladder) atomic system. The atomic states 0, 1 and 2 have 
transition frequencies U\ and U2- 
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FIG. 2: Time evolution of the probability of finding the system in state 2; P(t) = |&2(^)| 2 - The 
reservoir structure function is given by Eq. (|29j) with V =1 and = 1 in scaled units. The grid 
size for the discretised kernel was: 150x150 (solid), 100x100 (dashed), and 50x50 (dotted). In 
each case a range of ±30 for uj\ — luq and oo^ — loq was chosen. The result for a grid size of 150 x 150 
(solid curve) gives a reasonably accurate result. 
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FIG. 3: Time evolution of the probability of finding the system in state 2; P(t) = |&2(^)| 2 - The 
reservoir structure function is Eq. (|29[) with T =1 and: (a) O = 5.0; (b) 0, = 1.0; and (c) f2 = 0.3, 
in scaled units. The grid size for the discretised kernel was 150x150 chosen with a range of ±30 
for uj\ — u>o and uj^ — ujq in scaled units (as in Fig. EJ) . 
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FIG. 4: Time evolution of the probability of finding the system in state 2; P(t) = |&2(^)| 2 - The 
reservoir structure function is Eq. (|29j) with T =1 and £1 = 1.0. The two curves show the effect 
of changing from two separate reservoirs (solid line) to a single reservoir (dashed line). (Other 
parameters are as given in Fig. [2J The dashed line in this figure is identical to the solid line in Fig. 

El) 
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